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ABSTRACT 

Material-environment interactions are analyzed at 
microscopic scale to explain the lower than expected 
density observed by post-flight analysis of the char 
layer on the Stardust shield. Mass transfer, ablation 
(oxidation), and surface recession of fibrous material is 
simulated in 3D using a Monte-Carlo simulation tool. 
Ablation is found to occur either at the surface or in 
volume depending on Knudsen and Thiele number 
values. This study supports the idea of volume ablation 
followed by possible carbon fiber spallation that may 
explain post-flight analyses. 

INTRODUCTION 

Thermal protection systems (TPS) of planetary probes 
are exposed to extreme heat fluxes during 
atmospheric entry. The latest generation of TPS 
materials is fabricated by partial impregnation of 
carbon felt in phenolic resin. The material upper layer 
(typically 1 cm), exposed to temperatures higher than 
1200 K, is pyrolyzed into char. During pyrolysis, the 
phenolic resin loses up to 50% of its mass, shrinks 
and weakens. The pyrolyzed resin quickly oxidizes 
and may spall, revealing the structure of the felt. The 
felt can be described as a stack of carbon fibers 
randomly oriented (figure 1). The exposed carbon 
fibers progressively gasify by oxidation. This leads to 
mass loss of the carbon structure and wall recession 
called ablation. This work aims to improve the 
understanding of this latter process. 

Current models describe gasification of the char 
layer in terms of recession velocity (ablation velocity) 
of the overall surface. This description, based on 
models of previous generation TPS materials, has 
been shown to be correct for low porosity materials 
[Moyer, 1968]. However, the average surface may not 
recede uniformly in very porous materials (as in figure 
3), but the individual fibers may progressively vanish 
(as in figure 4) when oxidant penetration inside the 
porous material is fast enough compared to oxidant 
consumption by the fibers. In the second regime, 
ablation is no longer a surface phenomenon, but is 
volumetric. Seen from a surface point of view, two 
important consequences follow this volume ablation 
regime: (1) the effective reactivity of the material is 
significantly increased; (2) the material weakens in 
volume and the carbon fibers of the TPS material may 
be subject to mechanical erosion called spallation. In 
this context, the objective of this work is to understand 
when volume ablation may occur in order to 
determine whether spallation may have occurred 
during Stardust re-entry and to anticipate any 
spallation in future missions. 


As a first step in modeling volume ablation, this 
paper considers mass transfer and heterogeneous 
reaction coupling in random fibrous structure. An 
important point regarding mass transport during 
atmospheric entry is that wall pressure varies with 
atmospheric density and entry velocity. Hence, the 
molecular mean free path spans many orders of 
magnitude during entry leading to a wide span of the 
ablation regimes, successively, from Knudsen, to 
transition, and finally to bulk diffusion dominated 
mechanisms. 



Fig. 1 . Post-flight scanning electron microscopy 
micrograph of Stardust TPS surface [Stackpoole, 2008]. 


ORGANIZATION 

A phenomenological model for volume ablation is 
proposed. Three-dimensional (3D) numerical 
simulations of isothermal ablation in air of a carbon 
felt are performed and analyzed. To guide the 
analysis, a macroscopic model for volume ablation is 
derived analytically using volume averaging. The 
effective diffusivities needed to close the model are 
computed in bulk diffusion, transition and Knudsen 
regimes in porous media using direct numerical 
simulations. The model is used to provide 
understanding of the post-flight analysis of Stardust 
observations. 

MODEL 

Because the aim of this study is to improve the 
understanding of the char behavior, the focus is set 
on the material itself. Close to the wall, the velocity of 
the external airflow is null (no-slip condition); 
therefore mass transport can be considered diffusive. 
Knowledge of the value of the intensive variables at 
the wall (pressure, temperature, gas composition) is 
sufficient to simulate the material behavior. 

Carbon-based TPS is expected to ablate in air at 
a wall temperature T w around 3000K. We focus on 
the char layer that we assume to be composed of 
carbon fibers only, and that it ablates by oxidation. 
The layer is assumed to ablate along z, the normal 
direction, and is considered infinite in x and y 
directions (see Fig. 2). To start, the temperature of 
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the sample is assumed homogeneous and steady. In 
addition, blowing of pyrolysis gases is negligible. The 
structure of the carbon fibers is assumed 
homogeneous and isotropic. Finally, we assume 
heterogeneous oxidation reaction where carbon is 
consumed that leads to reduction of fiber radius. The 
motion of the interface between fiber and 
surrounding fluid can be interpreted as an advancing 
front with a normal velocity proportional to the 
oxidation rate fLachaud. 20081. 


c T p 

Wi r w 




a j TPS section b) Zoom on the top of the fiber stack 



- — 

c) Zoom on a section d) Elementary cell : capillary analogy 

Fig. 2. Multi -scale sketches of the computational cell 


The interface is represented by a surface function 
S(x,y,z,t) first order differentiable almost everywhere 
with constant value (zero) at the interface 
[Katardjiev, 1994]. The function S satisfies the 
differential equation: 

-gry 

— + v.VS = 0 (1) 

a 

where, 


v = Q-J-n (2) 

is the local surface velocity, with Q the solid molar 
volume, J the molar rate of oxidation, and 
n = VS / || VS || the normal pointing outwards from 

the surface. For a first order heterogeneous reaction, 
the local impinging molar flux density (on a fiber 
elementary surface) is given by: 


J = k f C (3) 

with k f the fiber intrinsic reactivity and C=C(x,y,z,t) 
the oxygen concentration. The local oxygen 
consumption is then coupled to the fiber surface 
recession and mass transfer into the porous medium. 
Mass transfer of oxygen in air is modeled by binary 


diffusion. Combining Fick’s law and mass 
conservation, oxygen transport writes: 

— + V-(-Z>VC) = 0 (4) 

dt 

The boundary conditions are: 

top (z=0): C = C W (Dirichlet), 

bottom (z=Ls): VC = 0 (Neumann), 
sides: periodic. 

For simple material architectures, this model can 
be solved analytically [Lachaud, 2008]. In the case of 
random fibrous media, direct numerical simulation or 
homogenization (averaging) methods are required. 
These solution methods are summarized in the two 
following sections. 

3D DIRECT NUMERICAL SIMULATION 

The fibers distribution at initial time is shown in the 
first pictures of Fig. 3 and 4. The computational 
domain is divided into the white part representing the 
fluid close to the wall and in the porous medium part, 
the fibers are represented by the black part. The 
porous medium is represented on a Cartesian grid in 
a cube of 100 3 voxels (3D pixels). Using a Monte 
Carlo algorithm, fibers with random positions and 
orientations are positioned in the cube until the 
porosity reaches the target value of 0.85. To 
reproduce realistically the material architecture, fibers 
are non-overlapping and their direction is azimuthally 
biased so that they are almost parallel to the 
topsurface (+/- 15°). 

The 3-D time-dependent solution to the 
reaction/diffusion models described above is obtained 
using an efficient numerical simulation code, named 
AMA. The code developed in previous work [Lachaud, 
2006], uses Monte-Carlo random walk whose 
ensemble averaged solution is represented by Eqs. 1 
and 4. AMA is a C ANSI implementation with four 
main features: (a) A 3-D image (graph) containing 
several phases (fluid and solid in the present case) is 
described by the discrete cubic voxels method on a 
Cartesian grid, (b) The moving fluid/solid interfaces 
are determined by a simplified marching cube 
approximation (Eq. 1). (c) Mass transfer by diffusion 
(Eq. 4) is simulated by a random walk using Maxwell- 
Boltzmann distribution for the velocity in the Knudsen 
and transition regimes, and a Brownian motion 
simulation technique in the bulk diffusion regime 
[Lachaud, 2006]. Brownian motion is a grid-free 
method that efficiently converges when simulating 
diffusion in a continuous fluid [Plapp, 2000]. (d) 
Fleterogeneous first-order reaction on the wall (Eq. 3) 
is simulated by using a sticking probability. Dirichlet 
boundary condition is handled using a buffer where 
the concentration of walker is kept constant. 

Two simulations carried out in bulk diffusion 
regime are shown in this paper (Fig. 3 and 4). The 
material behavior is found to depend on diffusion 
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velocity to oxidant consumption rate ratio inside the 
porous medium. When the reaction process is fast 
compared to mass transfer, the oxidant is mostly 
consumed at the surface and cannot seep into the 
porous medium; in this case ablation is a surface 
phenomenon (Fig. 3). In the converse case, the 
oxidant consumption being slow, the oxidant 
concentration becomes homogeneous in the porous 
medium and ablation occurs in volume (Fig. 4). This 
volume ablation obviously weakens the structure of 
the material, which is then likely to undergo spallation 
under the shear stress of the entry flow. At this point, 
an analytical study seems useful in order to find the 
relevant dimensionless number(s) for the involved 
coupling. 


fibers, the mean pore size writes [Transvalidou, 
1996]: 


i 46 * 

dp T 


( 5 ) 


For the material of interest in this study, d p is around 
50 fjm. For mean free paths lower than 1 pm 
Knudsen number (Kn=A/d p ) can be considered low 
and the regime continuous. In the validity domain of 
Maxwell-Boltzmann distribution, the mean free path 
of oxygen in air can be estimated using the following 
equation [Reid, 1987]: 


A = 95 - 10 ' 


-9 


io -r 

298 P 


( 6 ) 


ANALYTICAL SOLUTION (VOLUME AVERAGING) 
The behavior of the sample is analyzed at initial time, 
before ablation has significantly reduced the fiber 
diameter. 



Fig. 3. Diffusion regime: surface ablation. 


Fig.4. Reaction regime: volume ablation. 


In the literature, porous media are described in 
terms of mean pore size d p , porosity e, volume 
surface s, and tortuosity q. Assuming cylindrical 
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P (Pa) 

Fig.4. Validity domain of the continuous regime 
hypothesis. 


q z (z)-q 2 (z + dz) = q y (z) (10) 


-DV(C(z)-C(z + dz))S p = k f C(z)P f dz (11) 

where D is the local diffusion coefficient (D = D re f/r ] ), 
Pf is the local fiber perimeter, and Sp the local 
horizontal section of the pore. The integration of the 
ratio Sp/Pf on a representative horizontal section of 
the sample gives S/Pp e/s. In summary, C(z) is 
given by the following differential equation: 


5 2 C(z) 

dz 2 


4k, 


d p D 


C(z) = 0 


( 12 ) 


After integration, the oxygen concentration inside the 
porous medium is found to be: 


The temperature/pressure conditions for the validity 
of the continuous regime hypothesis are plotted in 
Fig. 4. Outside the validity domain, the mean free 
path is significantly lowered by gas collisions on the 
porous material walls. It results in a slower diffusive 
transport. Fick's law can be homogenized in all 
regimes. After homogenization in isotropic porous 
media, the law keeps the same form, but the 
effective diffusion coefficient writes [Tomadakis, 
1993]: 


C(z) = C 0 


cosh(0(z/L s -1) 
coshO 


where 


0 = 


Ls 



(13) 


(14) 


D „ = — D , 

e ff ^ ref 


(7) 


D ref is a reference diffusivity, corresponding to the 
longitudinal diffusivity into a capillary of diameter d p , 
which according to Bosanquet’s relation writes: 


— 1 — — =H = — (8) 

D ref D B D K Avd p 

with v the mean velocity of the molecules (molecular 

agitation) and X the mean free path. The tortuosity 
A], which is a geometric factor that characterizes the 
difference between a straight capillary and the actual 
tortuous medium as far as molecule trajectories are 
concerned, must be assessed numerically in the 
general case. 

Considering a first order heterogeneous reaction, 
the local impinging molar flux density (in an 
equivalent pore) is given by: 


q y (z) = k f C(z) 


(9) 


with C(z) the oxygen concentration at depth z (see 
Fig. 3-d). The mass balance in the elementary cell 
represented in figures 1-c and 1-d writes [Lachaud, 
2007]: 


is the Thiele number. The concentration gradient 
inside the porous medium is then a function of Thiele 
number, as shown on figure 5. When Thiele number 
is small, diffusion is much faster than reaction. In this 
case, the oxidant concentration inside the porous 
medium is homogeneous and equal to C 0 . Therefore, 
ablation is a volume phenomenon (Fig. 4). In the 
converse case (high Thiele number), diffusion is not 
fast enough to bring the oxidant inside the porous 
medium; only the top of the sample ablates. In this 
case, ablation is a surface phenomenon (Fig. 3). 



CM/C, 


In <t> 

3 

Fig.5. Oxidant concentration inside the porous medium as 
a function of Thiele number. 


In the design of a TPS, knowledge of the effective 
reactivity k e ff of the material is useful. k e ff is defined 
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as the reactivity of a homogeneous, dense and 
smooth material that would produce the same 
ablation rate as the porous medium: 

L s 

Jeff = Co^eff = J* C(z)k f sdz (15) 

z = 0 


walkers introduced randomly in the unit cell of the 
porous medium. The diffusion coefficient in direction 
j is then obtained using the following relation 
(Einstein, 1926): 


D -j = 



(17) 


The surface of the effective material is placed at 
abscissa z=0. The ratio of effective to fiber reactivity 
is obtained by integration of Eq. (15): 

k eff 4 eL v tanh^ „ 

t l = . , - +a-g) (i6) 

k f d p </> 



Fig. 6. Normalized effective reactivity as a function of the 
diffusion coefficient for kf=lm/s. 


The effective diffusivity inside the porous medium 

may be assessed for any combination ( v , /l ). 
However, equations (7) and (8) show that the only 
data that cannot be determined analytically is the 
tortuosity. 



Fig. 7. Path of a random walker (whose mean free path is 
5 micrometers) inside the porous medium of the study. 


As shown in Fig. 6, the material effective reactivity 
features two limit cases: surface ablation (Fig. 2) and 
volume ablation (Fig. 3). In the limit cases, the 
material behavior is no longer a function of the 
diffusion coefficient, while a strong dependency is 
observed in the transition regime. In the case of full 
volume ablation, the effective reactivity is found to be 
two orders of magnitude higher than the fiber 
intrinsic reactivity. 

EFFECTIVE DIFFUSIVITY (DNS) 

In order to apply this model to the analysis of actual 
problems, the knowledge of the effective diffusion 
coefficient inside the porous medium represented in 
figure 7 is required. This material being orthotropic, 
the diffusion coefficients in x-y and z directions are 
different. Effective Knudsen, transition, and Bulk 
diffusion coefficients in a given direction j are 
obtained through a computer simulation procedure 
using the mean square displacement in this 

direction, ^ 2 ) , after adequately large values of 

travel time, r , of a large population of random 


The tortuosity is only a function of X and of the 
porous medium architecture. As a consequence, it is 
convenient to assess the tortuosity as a function of 
Knudsen number (figure 8). The simulation tool 
presented in the first section is used again. In order 
to assure convergence, the simulations have been 
carried out with a total number of 10,000 walkers, 
and for a path one thousand times larger than the 
elementary cell side. The numerical procedure has 
been validated for a square array of parallel cylinders 
in comparison with analytical results [Perrins, 1979] 
in bulk diffusion regime and with numerical results in 
transition and Knudsen regime [Tomadakis, 1993; 
Vignoles, 2007]. Tomadakis and Sotirchos 
[Tomadakis, 1993] also studied random fiber 
structures but they were isotopic and with 
overlapping fibers, so they are not comparable to the 
anisotropic and non-overlapping fiber structure of 
this study. The extension of Bosanquet’s relation to 
ordinary porous media by defining a bulk tortuosity 
{ij B ) and a Knudsen tortuosity (r / K ) has been shown 
to provide a good interpolation of numerical data for 
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random fibrous structures [Tomadakis, 1993]. Under 
this extension, the effective tortuosity writes: 


' 1 + Kn 


( 18 ) 


The curves of figure 7 have been plotted using this 
approximation. A correct interpolation is obtained 
again for non-overlapping fibers. Relation (18) can 
then be used in the analytical ablation model. 


n 



Kn 

Fig. 8. Tortuosity of the random fibrous medium as a 
function of Knudsen number, for z and x-y directions 
(points: DNS; curves: analytical approximations). 


APPLICATION TO STARDUST ENTRY ANALYSIS 
At peak heating that corresponds to high ablation 
rate, pressure and temperature are around 10 4 Pa 
and 3000 K, respectively. In these conditions and in 
air, the mean free path is around 10 pm. The 
tortuosity in z direction (of interest for the ID 
analytical model) is then 1.23. Using Eq. (7), the 
diffusion coefficient is found to lie 

around 5-10 ~ 3 m 2 /s . In these conditions, the 
intrinsic reactivity of carbon fibers to air is known to 
be of the order of magnitude of 1 m/s [Drawin 1992, 
Lachaud 2007]. Figure 6 shows that the ablation 
regime is transitional and close to full volume 
ablation. The risk of spallation is high under the 
hypotheses of the model developed in this study. 
This explains the observation that the char density 
has been measured lower than expected during 
post-flight analyses [Stackpoole, 2008]. 

CONCLUSION 

A phenomenological model for the ablation of random 
fibrous architecture by oxidation has been proposed. 
The model includes oxidant mass transfer inside the 
porous medium, heterogeneous reaction, and fiber 
surface recession. Three-dimensional numerical 
simulations have been performed and have shown 
that ablation was controlled by the competition 
between reaction and diffusion and that it could be 
either a volume or a surface process. To help in the 


analysis, a simplified macroscopic model has been 
derived analytically using volume averaging. Thiele 
number has been found to be the key dimensionless 
number. When the reaction rate is fast compared to 
mass transfer rate (high Thiele number), the oxidant is 
mostly consumed at the surface and cannot seep 
through the porous char; in this case, ablation is a 
surface process. In the converse case (low Thiele 
number), the oxidant consumption rate being slow, 
the oxidant concentration becomes homogeneous 
inside the porous char and ablation occurs in volume. 
This volume ablation obviously weakens the structure 
of the material, which is then likely to undergo 
spallation under the shear stress of the entry flow. In 
order to feed the analytical model, the effective 
diffusivities in bulk, transition and Knudsen regime 
inside the porous medium have been computed by 
direct numerical simulation using a random walk 
algorithm. The tortuosity of the material as a function 
of Knudsen number is shown to follow a Bosanquet- 
like relation, from which the effective diffusivity can be 
computed. The model has finally been used to explain 
post-flight analysis of Stardust TPS ablation during its 
re-entry on Earth. Under the hypotheses of the model, 
ablation is found to be in volume at peak heating 
(corresponding to highe ablation rate). This explains 
the fact that the density of the material at the surface 
was lower that expected using surface ablation 
models. This implies the necessity of developing a full 
volume ablation model for a reliable design of the 
Crew Entry Vehicle TPS. This model should include 
pyrolysis gas advection, full chemistry, and heat 
transfer. 
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NOMENCLATURE 

Latin 

C oxygen concentration, mol/m 3 

d p mean pore diameter, m 

D diffusion coefficient, m 2 /s 
J molar ablation rate, mol/m 2 /s 

k reactivity, m/s 

Kn Knudsen number 

Ls sample length, or TPS thickness, m 
n vector normal to the surface 

P pressure, Pa 

P f local fiber perimeter, m 

q local molar flux density, mol/m 2 /s 

s volume surface, m 2 /m 3 

S p local horizontal section of a pore, m 2 

S(x,y,z,t) surface function 
T temperature, K 
v surface local velocity, m/s 
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v molecular agitation mean velocity, m/s 
x,y,z,t space (m) and time (s) coordinates 
Greek 

£ porosity 
r| tortuosity 

A mean free path, m 

mean square displacement in direction j, m 

p density, kg/m 3 

O Thiele number 

Q solid molar volume, m 3 /mol 
Subscripts 
B bulk 
e, eff effective 
f fiber 
K Knudsen 
ref reference 
w wall 
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